STOCHASTIC CELLULAR AUTOMATON FOR THE 
COAGULATION-FISSION-PROCESS 2A -» 3,4, 2A -> A 

HAYE HINRICHSEN 

Abstract. We introduce an efficient ceiluiar automaton for the coagulation- fission process 
with diffusion 2A — > 3A, 2A — > A in arbitrary dimensions. As the well-known Domany-Kinzel 
model, it is defined on a tilted hypercubic lattice and evolves by parallel updates. The model 
exhibits a non-equilibrium phase transition from an active into an absorbing phase and its 
critical properties are expected to be of the same type as in the pair contact process with 
diffusion. High-precision simulations on a parallel computer suggest that various quantities of 
interest do not show the expected power-law scaling, calling for new approaches to understand 
this unusual type of critical behavior. 



1. Introduction 

Non-equilibrium phase transitions from fluctuating phases into one or several absorbing states 
continue to fascinate many researchers because of their universal properties fl]-|3). The ultimate 
goal is to categorize such phase transitions into a finite number of universality classes. So 
far several universality classes have been identified, the most important ones being directed 
percolation (DP) the parity conserving (PC) class ||, voter-type transitions ||, and the 
generalized epidemic process [0]. Searching for other universality classes, the pair contact process 
with diffusion (PCPD) 2A — > 3A, 2A — > 0, also called annihilation-fission process, is one of the 
most promising candidates. In one spatial dimension the PCPD displays a novel type of critical 
behavior which does not seem to belong to any of the known universality classes. In contrast to 
DP and PC transitions, where individual particles can generate offspring, the PCPD is a binary 
spreading process, i.e., two particles have to meet in order to create a new one. 

The possibility of a non-trivial critical behavior in binary spreading processes was already 
pointed out by Grassberger in 1982 ||]. Nethertheless it took another 15 years before the PCPD 
was investigated systematically for the first time by Howard and Tauber ||. Their work was 
motivated by the observation that the annihilation process 2A — * produces anticorrelations 
among the particles (emerging as 'imaginary' noise in the Langevin equation), while the fission 
process 2 A — > 3 A generates positive correlations ('real' noise). Combining the two reactions 
they were able to show that the interpolation between 'real' and 'imaginary' noise leads to a 
non-equilibrium phase transition. However, the corresponding bosonic field theory turned out 
to be non-renormalizable. This is partly due to the fact that in the bosonic formulation the 
density of particles diverges in the active phase so that perturbative methods cannot be applied. 
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Carlon et al. |T(| were the first to consider a fermionic lattice model of the PCPD, in which 
the local density of particles is bound by an exclusion principle. Their paper triggered a series of 
numerical studies |TT|-[r^| and released a still ongoing debate concerning the critical behavior at 
the transition. Even today, two years later, it is not yet clear to what extent the PCPD represents 
an independent genuine universality class. On the one hand, a fermionic field theory for the 
PCPD does not yet exist. On the other hand, density matrix renormalization techniques and 
Monte Carlo simulations turned out to be unreliable because of strong deviations from power law 
scaling. Moreover, in some studies the estimated critical exponents seem to depend on the choice 
of the diffusion rate, leading Odor to the conjecture that the transition may be characterized by 
two independent universality classes for weak and strong diffusion Reported estimates for 
the critical exponents are scattered over a wide range, depending on the numerical effort and 
the value of the diffusion rate. Nevertheless there is a consensus that the PCPD does exhibit 
a novel type of non-equilibrium phase transition with a very characteristic critical behavior. 
The transition is universal in the sense that the same characteristic critical properties can be 
observed in various other binary spreading processes such as in multi-component models [ [17|| , 



parity-conserving variants | |15| , and especially the coagulation-fission process [14,18] 
(1) 2A^3A, 2A^A, 

on which we will focus in the present work (see Fig. [I]). By analyzing this process we will 
implicitely assume that the conclusions can also be applied to the PCPD and other binary 
spreading processes, at least on a qualitative level. 

What is the reason for the unusually strong deviations from power-law scaling in the 1+1- 
dimensional PCPD? Initially it was argued that the deviations may be due to logarithmic correc- 
tions at the upper critical dimension. However, Odor showed that the upper critical dimension 



of the process is d c = 2 [16|. For this reason ordinary logarithmic corrections at the upper 



critical dimension can be ruled out unless the system has several critical dimensions. 

In order to explain the observed phenomena phenomenologically, it was suggested that the 
PCPD can be interpreted as a DP process coupled with an annihilation process in a feedback 
loop Numerical simulations confirmed that such a cyclically coupled model displays in- 

deed a similar characteristic critical behavior at the transition. Based on these ideas Noh and 
Park |l^] proposed that the transition in the PCPD might be characterized by continuously 
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varying exponents depending on the value of the diffusion rate. As an explanation they sug- 
gested that the PCPD may be viewed as a DP process subjected to a marginal perturbation 
caused by long-term memory effects of the dynamics. 

Even more recently Park and Kim presented numerical evidence that the PCPD is char- 
acterized by a well-defined set of critical exponents provided that the rates for diffusion and 
annihilation (or coagulation) are tuned in such a way that the process without branching be- 
comes exactly solvable |20|. Very recent studies even suggest that all PCPD variants eventually 
cross over to the same asymptotic critical behavior with a unique set of critical exponents 
meaning that the PCPD would indeed represent an independent universality class. 

So far all explanations anticipate that the critical PCPD displays asymptotic power laws. 
However, in the present study we find that there is no power-law scaling at all, at least within 
the numerically accessible range. To this end we introduce a very simple and efficient stochastic 
cellular automaton model for the coagulation- fission process (|l|). Since the model evolves by 
synchronous updates, it can be implemented easily on a parallel computer. This allows us 
to perform high-precision simulations exceeding the number of time steps used in previous 
studies p3|. Our numerical results suggest that even after one million time steps the scaling 
regime is not yet reached. Inspired by this observation we discuss the possibility that the process 
may display an extremely slow crossover to directed percolation. 



2. Definition of the cellular automaton 

The cellular automaton for the 1+1-dimensional coagulation-fission process investigated in the 
present work is defined as follows. As the well-known Domany-Kinzel model [23], it uses a tilted 



square lattice, where the sites are enumerated horizontally by a position index i and vertically 
by a temporal index t. Since the lattice is tilted the position indices i at even (odd) times t are 
also even (odd). The sites of the lattice can be either active (sj(i) = 1) or inactive (si(t) = 0). 
Each t labels a horizontal row, specifiying the state of the model at time t. 

The model evolves by parallel updates as follows. Assume that the state at time t is given. 
In order to construct the new state at time t + 1, we first clear all sites of the next horizontal 
row. Then a loop is executed which runs over all active sites j of the previous configuration at 
time t. For each of them we first check whether its nearest neighbors Sj-2(t) and Sj+2(t) are 
occupied. Depending on their state the following local updates are carried out: 

(a) If both neighbors are inactive, the particle performs a coagulating random walk by acti- 
vating either site s,_ i(t + 1) or Sj+i(t + 1) in the new configuration with equal probability. 

(b) Otherwise, if at least one of the neighbors is active, another random number z € [0, 1] 
is generated. If z < p the two connected sites are activated by setting Sj-±(t + 1) = 
Sj + i(i + 1) = 1, whereas for z > p either site s,_ i(t + 1) or Sj+i(i + 1) is activated with 
equal probability in the same way as in (a). 

It is easy to see that these dynamic rules resemble the process 2 A — > 3 A, 2 A — > A combined with 
diffusion. Unlike in random-sequential models, it is not necessary to select particles randomly 
so that less random numbers have to be generated. A minimal pseudo code for the update 
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Figure 2. Dynamic rules of the cellular automaton in 1+1 dimensions, (a) 
Particles without a nearest neighbor perform a coagulating random walk by ac- 
tivating one of the adjacent sites in the new configuration with equal probablity. 
(b) If a particle has at least one nearest neighbor (indicated by grey color), both 
following sites are activated with probability p. Otherwise, the particle performs 
a coagulating random walk as in (a). 



algorithm is given in the appendix. In order to simulate the cellular automaton efficiently, 
further optimizations are necessary. Instead of using a static array describing a finite lattice 
with periodic boundary conditions, it is useful to store the positions of active sites at time t 
individually in a dynamically generated list. Using such a list it is no longer necessary to perform 
a loop over all sites and to check whether they are active, instead the loop runs exclusively over 
the active sites, speeding up the simulation especially when the particle density is low. List-based 
algorithms are also advantageous in so-called dynamical simulations starting with a localized 
active seed p3 . Here the system size is practically infinite (limited only by the range of integer 
numbers so that finite-size effects can be excluded. 

Moreover, the update procedure can be implemented easily on a parallel computer. To this 
end we note that the new state at time t + 1 is initially set to zero and then constructed by an 
unconditional activation, i.e., we assign the value 1 to certain sites irrespective of their state. 
For this reason the total update can be understood as a Boolean 'OR' operation of all local 
updates within the loop. This means the local updates commute and thus can be processed 
in parallel blocks. Therefore, it is straight forward to implement the cellular automaton on a 
parallel computer with a ring geometry. Finally we note that the model can easily be generalized 
to higher dimensions in the same way as the Domany-Kinzel model, although the present work 
will be restricted to the 1+1-dimensional case. 



3. Numerical results 

3.1. Homogeneous initial conditions. Fig. || shows the numerical results for the decay of 
the particle density close to the transition starting with a fully occupied lattice. In systems 
with power-law scaling the density is expected to decrease as p(t) ~ t~ s , where 5 = (3/vu. The 
parallelized implementation allows us to simulate very large systems with 2 21 «2x 10 6 sites so 
that finite-size effects can be excluded. The simulation time of 2.5 x 10 6 updates exceeds the 
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Figure 3. Coagulation-fission process in 1+1 dimensions: Decay of the particle 
density close to criticality starting with homogeneous initial conditions. The left 
graph shows the raw data of seven independent simulations for p running from 
0.29992 (bottom) to 0.30004 (top) in steps of 0.00002. The right graphs shows 
the same data multiplied by t 0,246 . 

range of previous numerical works. Because of the unusually large system size it is sufficient to 
average the density of particles over 32 independent runs. 

In order to estimate the critical threshold, we measured the decay for different values of p. 
Plotting the raw data in a double logarithmic plot (left panel of Fig. ||), there seems to be a 
reasonable straight line with a slope —5 ~ —0.246(10) over almost four decades. However, if the 
measured values are multiplied by t - 246 (right panel), it can be seen that the lines are actually 
curved. Note that the data sets for different p are statistically independent so that the observed 
trend of the curvature cannot be attributed to correlated fluctuation effects. 

The figure illustrates how difficult it is to determine the critical point. If the curves eventually 
veer down, we assume that the process is subcritical. Using this criterion we find that p c > 
0.29998. The exact critical point, which is presumably somewhere between 0.3000 and 0.3002, 
cannot be determined precisely since all data sets for p > 0.29998 show a non-vanishing positive 
curvature over at least three decades. Thus we have strong reasons to believe there no power-law 
scaling, at least in over the first six decades in time. In less extensive simulations one would 
have been tempted to interpret the broad plateau between 10 2 and 10 4 time steps as the onset of 
a power-law and to lower the percolation probability in order to get a straight line. The present 
numerical results show that this is not possible. 

3.2. Seed simulations. Studying phase transitions into absorbing states it is customary to 
perform simulations starting with a localized seed [Q. The quantities of interest are the survival 
probability P(t), the number of particles N(t) averaged over all runs, and the mean square 
spreading from the origin R 2 (t) averaged over the surviving runs. In the standard scaling theory 
of phase transitions into absorbing states these quantities are expected to obey asymptotic power 
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Figure 4. Seed simulations. Left: Survival probablity as a function of time 
averaged over 6 x 10 6 realizations for p running from 0.29991 (bottom) to 0.30012 
(top) in steps of 0.00003. Right: The same data multiplied by t ' 12 . 




t [MCS] t [MCS] 

Figure 5. Seed simulations. Left: Number of particles (upper curves) and 
number of particles with at least one active nearest neighbors (lower curves) as 
a function of time. Right: Mean square spreading from the origin averaged over 
surviving runs. 



laws of the form 

(2) P(t)~t~ 6 ', N(t)~f>, R 2 {t)~t 2/z . 

In the coagulation-fission process seed simulations are conceptually difficult. Since the absorbing 
state corresponds to a single diffusing particle, one has to start with a localized pair of particles. 
The survival probability is then most naturally defined as the probability that the system has 
not yet reached the absorbing state, i.e., the process is called 'surviving' as long as there are at 
least two particles in the system. However, this definition is problematic in so far as it is no 
longer possible to divide the system into active sites and absorbing domains [25]. Irrespective 
of these difficulties, the survival probability P(t) displays a similar behavior as the decay of the 
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exponent 


t m 10 ' 2 


t 10 b 


DP-value 


5 


0.25(1) 


0.22(1) 


0.159 


5' 


0.11(1) 


0.13(1) 


0.159 


V 


0.07(2) 


0.22(3) 


0.314 


z 


2.00(4) 


1.78(5) 


1.581 



Table 1 . Effective exponents moving towards DP values as time proceeds. 



density in the case of homogeneous initial conditions. Plotting the raw data in a log- log plot 
one might expect a straight line with an approximate slope of —5' = —0.12(1) (see left panel 
of However, multiplying the data points by i ' 12 (right panel of Fig. 0) there is again a clear 
curvature over more than four decades. Thus, we can exclude power-law scaling, at least in the 
numerically accessible range. 

Similar results are obtained for the other quantities. The left panel of Fig. || shows the number 
of particles N(t) averaged over all runs. In addition, the figure shows the average number of 
particles with at least one active nearest neighbor A^i), i.e., the number of particles which can 
generate offspring. As previously observed in other studies of binary spreading processes in 1+1 
dimensions, both quantities behave similarly and their quotient seems to tend to a constant. 
However, also in this case it is obvious that they do not follow a power-law. Likewise, the mean 
square spreading from the origin R 2 (t) does not display power-law scaling. As shown in the 
right panel of Fig. ||, after an initial transient there is horizontal plateau, where R 2 (t) increases 
linearly, indicating an intermediate diffusive behavior. After 10 3 time steps, however, the curves 
veer up. This means that the effective exponent z{t) decreases with time so that the process 
spreads superdiffusively. 



4. Does the PCPD belong to the DP class? 

What is the origin for the unusual critical behavior observed in the coagulation-fission process. 
Why is it impossible to obtain clean power laws? One possibility, which has not been discussed 
so far, would be that the critical behavior of the PCPD and related models crosses over to DP 
after very long time. Several arguments support this speculative point of view: 

1. Unlike unary spreading processes, which generate offspring everywhere along the particle 
trajectories, the PCPD and related models are binary spreading processes, where solitary 



particles can move diffusively over long distances. As pointed out in Ref. [13|, the dynamics 
switches between two modes, namely, dense regions in space-time, where branching occurs 
frequently, and sparse regions, where solitary particles diffuse over long distances. The 
apparent contradiction that single particles diffuse (with a dynamic exponent z = 2) while 
the whole process is super diffusive (z < 2) can only be resolved if the effective diffusion 
rate itself scales to zero under rescaling in the renormalization group sense. Therefore, the 
critical PCPD regarded on extremely large scales should look like a pair contact process 
without diffusion, which is known to belong to the DP universality class [26|. 
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2. As shown in Table El all effective exponents move in the direction of the DP values as time 
proceeds. The same applies to the density exponent [3 (not reported here). However, the 
accuracy of the numerical results does not permit a reliable extrapolation. 

3. As shown in Fig. ||, the effective dynamic exponent z deviates only slightly from 2. Between 
10 and 100 time steps there is virtually no deviation. At 10 4 time steps the deviation is 
less then 8% and after 10 6 time steps it is still less than 12%. The small deviation makes it 
plausible why the diffusion constant changes only gradually under renormalization, leading 
to an extremely slow crossover. 

4. In the inactive phase a binary spreading process in 1+1 dimensions displays an algebraic 
decay p(t) ~ i -1 / 2 , whereas in DP the density is known to decay exponentially. This 
apparent contradiction can be resolved by observing that the decay of the particle density 
in a slightly subcritical PCPD can be divided into three different temporal regimes. In the 
first regime (region I in Fig. ^) spreading and diffusion are still strong enough to sustain 
each other so that the particle density decays very slowly. After a characteristic time, 
however, the density decays quickly (region II) until it crosses over to an algebraic decay 
£-1/2 ( re gi on Tjj) Remarkably, the amplitude of this asymptotic power law does not depend 
on p. Consequently, the sudden jump in region II becomes more and more pronounced as 
the critical threshold is approached and may tend to a quasi-exponential decay in the limit 

P^Pc- 

In spite of these arguments the hypothesis of an eventual crossover to DP faces several problems. 
If it were valid in higher dimensions, it would be in conflict with Odors observation |l(| that 
the upper critical dimension of the PCPD is smaller than 4. Moreover, it is not clear how the 
time reversal symmetry of DP is recovered in the asymptotic limit. In addition, the implicit 
assumption that all variants of binary spreading processes behave essentially in the same way 
has to be questioned. Finally, other high-precision simulations are reported that seem to confirm 
power-law behavior f2l"| , |22|| . 
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5. Summary 



In the present work we have introduced an efficient cellular automaton model for the coagulation- 
fission process with diffusion. Performing extensive simulations on a parallel computer, we find 
that the model does not exhibit power-law scaling in the numerically accessible range, instead 
the effective exponents (local slopes) show a continuous drift. As in various other studies, we 
may have been tempted to produce almost straight lines over the last few decades by adjusting 
the critical parameter, but the overall curvature of the data sets suggests that in this case the 
system would be slightly subcritical. For this reason we conjecture that the drift of the effective 
exponents will continue far beyond the numerically accessible range. We expect the same to be 
true for other binary spreading processes such as the PCPD. 

Observing that all effective exponents move towards the corresponding DP exponents (al- 
though they are still quite far away), one has to consider the possiblity that the system may 
cross over to DP after very long time. As a possible explanation we suggest that the diffusion 
rate may scale to zero in the renormalization group sense, driving the system eventually to a 
binary spreading process without diffusion, which is expected to display the critical behavior 
of DP. To this end it would be helpful to perform a real-space renormalization group study in 
order to see how the parameters, especially the diffusion rate, change under rescaling p8[. 



To summarize, there is a confusing variety of opinions concerning the critical behavior of 
binary spreading processes. Presently, the main viewpoints are: 



- Binary spreading processes belong to a single universality class with a well-defined set of 



critical exponents [20|. 



Numerical results indicate that there are two independent universality classes for low and 



high values of the diffusion rate [10, 12|. 



The PCPD may be interpreted as a cyclically coupled DP and annihilation process [13]. 



The PCPD can be regarded as a marginally perturbed DP process with continuously vary- 
ing exponents fl9[ . 

Binary spreading processes may cross over to DP after very long time (present work). 
Finally, one cannot exclude the possibility that the PCPD and related models exhibit a 



weakly discontinuous transition [27] 



Further investigations are needed to find out which of these explanations is the correct one. 
Finally we note that the same questions arise in the context of a recently introduced triplet 
process with diffusion 3A — > 4A, 3A — > [29]. As in the case of the 1+1-dimensional PCPD, it 
may be possible that this process crosses over to DP as well after very long time. 
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Appendix A. Minimal pseudo code 

A minimal pseudo code (disregarding the boundary conditions) for the update procedure of the 
cellular automaton can be written as: 

integer L,T; // system size and simulation time 
integer s[L,T]; // the lattice 
float p; // the branching probability- 
integer j ; // site index 
integer start = t mod 2; // select even/odd sublattice 

for j from 1 to L-l do s[j,t+l]=0; // clear new state at time t+1 

for j from start to L-l in steps of 2 do // loop over all active sites 
if s [j ,t] ==1 then 

if random<0.5 then s [j-1 ,t+l] =1 ; else s [j+1 ,t+l] =1 ; 
if s[j-2,t]+s[j+2,t]>0 and random<p then s [j-1 ,t+l] =s [j+1 ,t+l] =1 ; 
end (if) 
end (for) 



stochastic cellular automaton for the coagulation-fission-process 2a -* 3a, 2a -> all 
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